function ImageBinary=PlotOverlap_angle_L145(L1,L4,L5,epoch,FN_i,CameraAngleCapability,DU,VU,TU)

options=odeset('RelTol',1e-13,'AbsTol',1e-13);
[xs,ys,zs] = ellipsoid(-0,0,0,696340,696340,696340,500);


%% Compute L4

S_i = L4.S_i;

temp_time_month = month(L4.epoch);

ind = find(temp_time_month==epoch);

% plot3(S(1,:),S(2,:),S(3,:))
[row,col,~] = size(FN_i);
visibility_i=zeros(row,col);
parfor ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if angle_temp<CameraAngleCapability
                visibility_i(ii,jj) = angle_temp;
            else
                visibility_i(ii,jj) = nan;
            end
        end
    end
end


ImageBinary_L4 = visibility_i;

% figure;
% hold on;
% mesh(ImageBinary_L4*100,'EdgeColor','none','FaceColor','flat');
% xlabel('Longitude')
% ylabel('Latitude')
% titlestring = sprintf('L4 Month = %d',epoch);
% title(titlestring)
% set(gca,'fontsize',16)
% xticks(linspace(0,row,5));
% xticklabels({'180W','90W','0','90E','180E'});
% yticks(linspace(0,col,5))
% yticklabels({'90S','45S','0','45N','90N'});
% xlim([0,row]);
% ylim([0,col]);

%% Compute L5

S_i = L5.S_i;

temp_time_month = month(L5.epoch);

ind = find(temp_time_month==epoch);

% plot3(S(1,:),S(2,:),S(3,:))
[row,col,~] = size(FN_i);
visibility_i=zeros(row,col);
parfor ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)))
            if angle_temp<CameraAngleCapability
                visibility_i(ii,jj) = angle_temp;
            else
                visibility_i(ii,jj) = nan;
            end
        end
    end
end

S_i = L5.S_i;

temp_time_month = month(L5.epoch);

ind = find(temp_time_month==epoch);

[row,col,~] = size(FN_i);
visibility_i_L5_face=zeros(row,col);
parfor ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)))
            if and((89<angle_temp),(angle_temp<91))
                visibility_i_L5_face(ii,jj) = 100;
            else
                visibility_i_L5_face(ii,jj) = nan;
            end
        end
    end
end

ImageBinary_L5_face = visibility_i_L5_face;


ImageBinary_L5 = visibility_i;

% figure;
% hold on;
% mesh((ImageBinary_L5+ImageBinary_L5_face)*100,'EdgeColor','none','FaceColor','flat');
% xlabel('Longitude')
% ylabel('Latitude')
% title('L5')
% set(gca,'fontsize',16)
% xticks(linspace(0,row,5));
% xticklabels({'180W','90W','0','90E','180E'});
% yticks(linspace(0,col,5))
% yticklabels({'90S','45S','0','45N','90N'});
% xlim([0,row]);
% ylim([0,col])

%% Compute L1

S_i = L1.S_i;

temp_time_month = month(L1.epoch);

ind = find(temp_time_month==epoch);

[row,col,~] = size(FN_i);
visibility_i=zeros(row,col);
for ii = 1:row
    for jj = 1:col
        normvec = [FN_i(ii,jj,1);FN_i(ii,jj,2);FN_i(ii,jj,3)];
        for k = ind(1)
            posvec = S_i(1:3,k);
            angle_temp = acosd(dot(normvec,posvec)/(norm(normvec)*norm(posvec)));
            if angle_temp<CameraAngleCapability
                visibility_i(ii,jj) = angle_temp;
            else
                visibility_i(ii,jj) = nan;
            end            
            if round(angle_temp)==0
                earthLocation=[ii,jj];
            end
        end
    end
end


ImageBinary_L1 = visibility_i;
% figure;
% hold on;
% mesh(ImageBinary_L1*100,'EdgeColor','none','FaceColor','flat');
% xlabel('Longitude')
% ylabel('Latitude')
% title('L1')
% set(gca,'fontsize',16)
% xticks(linspace(0,row,5));
% xticklabels({'180W','90W','0','90E','180E'});
% yticks(linspace(0,col,5))
% yticklabels({'90S','45S','0','45N','90N'});
% xlim([0,row]);
% ylim([0,col]);
% colormap("gray");
%% L4+L5

% ImageBinary_L145 = ImageBinary_L1+ImageBinary_L5+ImageBinary_L4+ImageBinary_L5_face;

figure;
hold on
contour(ImageBinary_L4,'LevelList',10:20:70,'ShowText','on');
contour(ImageBinary_L1,'LevelList',10:20:70,'ShowText','on');
contour(ImageBinary_L5,'LevelList',10:20:70,'ShowText','on');
% graycolor = gray(10);
colormap("turbo")
% mesh(ImageBinary_L5_face,'EdgeColor','interp','FaceColor','interp');
xlabel('Longitude')
ylabel('Latitude')
titlestring = sprintf('L4 Month = %d',epoch);
% title(titlestring)
set(gca,'fontsize',15)
xticks(linspace(0,row,5));
xticklabels({'180W','90W','0','90E','180E'});
yticks(linspace(0,col,5))
yticklabels({'90S','45S','0','45N','90N'});
xlim([0,row]);
ylim([0,col])
plot3(earthLocation(2),earthLocation(1),300,'ro','MarkerSize',7,'MarkerFaceColor','r')



end